################################################################################
# This script was created to recreate microclimate based on a linear model
# This script will read in some static ASCII data (coastdist, EVI)
# And loop through Tmax and Solar data from AWAP

################################################################################
#load necessary libraries

#list the libraries needed
necessary=c("adehabitat","SDMTools","rgdal","sp")

#check if library is installed
installed = necessary %in% installed.packages()

#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)

#load the libraries
for (lib in necessary) library(lib,character.only=T)

################################################################################

# Define Tmax Directory
tmax.dir = "/home1/99/jc152199/MicroclimateStatisticalDownscale/SpatialTMax/AWAP/Tmax/2007/"

# Define Solar Directory

solar.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/SOLAR/'
solar.join = data.frame(juldate=c(dir(solar.dir)), month=c(1:12))

# Define the out directory
 
out.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/MicroclimateASCII/Max/'

# Define the static directory

static.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/STATIC/'

# Create a vector that is coefficients from the linear model

coefs = as.numeric(c(10.2443755728, 1.0147051738,-0.0730804885, 0.0012472567,-0.0002502741,-0.0007090372))
coef.names =c('Intercept','AWAP_max','fpcmean','fpcvar','solar','roaddist')
names(coefs)=coef.names

# Read in the static ASCII data (coastdist and EVI)

fpcmean.asc = read.asc(paste(static.dir,'fpcmeanwtplusbuffer250m.asc',sep=""))
fpcvar.asc = read.asc(paste(static.dir,'fpcvarwtplusbuffer250m.asc',sep=""))
#coastdist.asc = read.asc(paste(static.dir,'coastdistwtplusbuffer250m.asc',sep=""))
roaddist.asc = read.asc(paste(static.dir,'disttoroadwtplusbuffer250m.asc',sep=""))

# Define a blank ASCII (all cells zero) that is the same extent and cell size as AWAP data need to be 
base.asc = roaddist.asc 
base.asc = base.asc * 0

## This is the base.pos command as given to me by JJV.  It only identifies row/column combo's which store NA as their value 
### This command was used to generate my spatial surfaces of Tmax
base.pos = as.data.frame(which(is.finite(base.asc), arr.ind = T)) 


# The follow command creates a two-column data frame whose number of rows is equal to the number of cells in the file base.asc
# Row/column position is retained

base.pos = expand.grid(c(1:nrow(base.asc)),c(1:ncol(base.asc)))
base.pos = data.frame(X=base.pos[,1],Y=base.pos[,2])

# Add east and north to this data frame of row/column positions

base.pos$east = getXYcoords(base.asc)$x[base.pos$X]
base.pos$north = getXYcoords(base.asc)$y[base.pos$Y]

# Convert east-north to lat-long
### EPSG 32755 WGS1984 Lat Long
### EGSG 4326 WGS1984 UTM

tout = as.data.frame(spTransform(SpatialPoints(cbind(base.pos[3:4]),proj4string=CRS("+init=epsg:32755")), CRS("+init=epsg:4326")))
base.pos$lon = tout[,1]
base.pos$lat = tout[,2]

# Create a matrix with the same rows and columns as .asc file
# The value of each cell will be the coefficient value of that independent variable

fpcmean.mat = matrix(coefs[3],nrow(base.asc),ncol(base.asc))
fpcvar.mat = matrix(coefs[4],nrow(base.asc),ncol(base.asc))
roaddist.mat = matrix(coefs[6],nrow(base.asc),ncol(base.asc))

# Multiply the ASCII files by the matrices and create a 'static' ASCII layer that will be added into the following loop

static.asc = (roaddist.asc * roaddist.mat) + (fpcmean.asc * fpcmean.mat) + (fpcvar.asc * fpcvar.mat)

# Define the year

year = 2007

  #doy = 0 #track the number of days in the year
  for (month in 1:12) {
                   
    cat(month,'\n')                      
    
    #extract solar data for the month in question
    t.join = solar.join[which(solar.join$month==month),1] 
    solar.asc = read.asc(paste(solar.dir,solar.join[which(solar.join$month==month),1],sep=""))
    solar.mat = matrix(coefs[5],nrow(solar.asc), ncol(solar.asc))
    solar.coef.asc = solar.asc * solar.mat
    
    #dom = 0 #track the number of days in the month
    for (day in 1:31) {
      
      #calculate the date
      tdate = format(as.Date(paste(year,month,day,sep="-"),"%Y-%m-%d"),"%Y%m%d")
      cat(tdate,'\n')
      
      #if the date exists as a file... start the summary
      tfile = paste(tmax.dir,"tmax.",tdate,tdate,".grid.gz",sep="")
      if (file.exists(tfile)){  cat(tfile,'\n')
      
        # import the asc file  and define the multiplication matrix
        awap.asc = read.asc.gz(tfile)
        awap.mat = matrix(coefs[2], nrow(base.asc), ncol(base.asc))
        
        # Extract data from AWAP ASCII, storing maxtemp against lat-long
        
        awap.data = extract.data(cbind(base.pos$lat,base.pos$lon),awap.asc)
        awap.pos = cbind(base.pos$lon, base.pos$lat, awap.data)
                                            
        # Define a blank ASCII in the loop, populate it with values from AWAP extract
                                     
        t.asc = base.asc 
        t.asc[cbind(base.pos$row,base.pos$col)] = awap.pos[,3]
        base.asc[cbind(base.pos$row,base.pos$col)]=t.data[,5]
        
        # Multiply AWAP ASCII by coef matrix
        
        awap.coef.asc = t.asc * awap.mat
        
        # Define microclimate as an ASCII sum of static/solar/awap/intercept
        
        final.asc = static.asc + solar.coef.asc + awap.coef.asc + coefs[1]
        cat(head(final.asc),'\n')
        
        # Write out an ASCII file of max temp
        
        write.asc.gz(final.asc, paste(out.dir,tdate,"_Tmax",sep="")) 
        
          }
        }
      }


